function ivol = blsImpVolCall( S, K, r, T, price, q, tolerance )

N_iterations = 10;

%S, K and price are column vectors. r, T and q are scalars
xn = 0.2*ones(length(S),1);
xn(isnan(price)) = NaN;

for i = 1:N_iterations
    temp1 = ( log(S./K)+r-q )./( xn.*sqrt(T) );
    d1 = temp1 + xn.*sqrt(T)/2;
    d2 = temp1 - xn.*sqrt(T)/2;

    phi1 = exp( -d1.^2 / 2 ) ./ sqrt(2*pi);
    
    f = S.*exp(-q*T).*normcdf(d1) - K.*exp(-r.*T).*normcdf(d2) - price;
    ind = find( abs(f) > tolerance );
    if isempty(ind)
        ivol = xn;
        return;
    end
     xn(ind) = xn(ind) - f(ind) ./ (sqrt(T(ind)).* phi1(ind).*S(ind).*exp(-q*T(ind)) );
     
    if i==N_iterations
        ivol = xn;
        ivol(ind) = NaN;
        return
    end
end
